CB624

Pressure Trajectories

Show code
library(GeoPressureR)
library(leaflet)
library(leaflet.extras)
library(raster)
library(dplyr)
library(ggplot2)
library(kableExtra)
library(plotly)
library(GeoLocTools)
setupGeolocation()
knitr::opts_chunk$set(echo = FALSE)
load(paste0("../data/3_static/", params$gdl_id, "_static_prob.Rdata"))
load(paste0("../data/4_basic_graph/", params$gdl_id, "_basic_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_wind_graph.Rdata"))
load(paste0("../data/5_wind_graph/", params$gdl_id, "_grl.Rdata"))
col <- rep(RColorBrewer::brewer.pal(8, "Dark2"), times = ceiling(max(pam$sta$sta_id) / 8))

Shortest path and simulated path

The large circles indicates the shortest path (overall most likely trajectory) estimated by the graph approach. The size is proportional to the duration of stay. The small dots and grey lines represents 10 possible trajectories of the bird according to the model.

Click on the full-screen mode button on the top-left of the map to see more details.

Show code
sta_duration <- unlist(lapply(static_prob_marginal, function(x) {
  as.numeric(difftime(metadata(x)$temporal_extent[2], metadata(x)$temporal_extent[1], units = "days"))
}))
pal <- colorFactor(col, as.factor(seq_len(length(col))))
m <- leaflet(width = "100%") %>%
  addProviderTiles(providers$Stamen.TerrainBackground) %>%
  addFullscreenControl()

for (i in seq_len(nrow(path_sim$lon))) {
  m <- m %>%
    addPolylines(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = 0.5, weight = 1, color = "#808080") %>%
    addCircles(lng = path_sim$lon[i, ], lat = path_sim$lat[i, ], opacity = .7, weight = 1, color = pal(factor(shortest_path$sta_id, levels = pam$sta$sta_id)))
}
m %>%
  addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
  addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = pal(factor(shortest_path$sta_id, levels = pam$sta$sta_id)), weight = sta_duration^(0.3) * 10, popup = paste0("sta_id=", shortest_path$sta_id))

Marginal probability map

The marginal probability map estimate the overall probability of position at each stationary period regardless of the trajectory taken by the bird. It is the most useful quantification of the uncertainty of the position of the bird.

Show code
li_s <- list()
l <- leaflet(width = "100%") %>%
  addProviderTiles(providers$Stamen.TerrainBackground) %>%
  addFullscreenControl()
for (i_r in seq_len(length(static_prob_marginal))) {
  i_s <- metadata(static_prob_marginal[[i_r]])$sta_id
  info <- metadata(static_prob_marginal[[i_r]])$temporal_extent
  info_str <- paste0(i_s, " | ", info[1], "->", info[2])
  li_s <- append(li_s, info_str)
  l <- l %>%
    addRasterImage(static_prob_marginal[[i_r]], colors = "OrRd", opacity = 0.8, group = info_str)
}
l  %>%
  addLayersControl(
    overlayGroups = li_s,
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  hideGroup(tail(li_s, length(li_s) - 1))

Altitude

Altitudes are computed based on pressure measurement of the geolocation, corrected based on the assumed location of the shortest path. This correction accounts for the natural variation of pressure as estimated by ERA-5.

Show code
p <- ggplot() +
  # geom_line(data = pam$pressure, aes(x = date, y = obs), colour = "grey") +
  geom_line(data = do.call("rbind", shortest_path_timeserie), aes(x = date, y = altitude)) +
  geom_line(data = do.call("rbind", shortest_path_timeserie) %>% filter(sta_id > 0), aes(x = date, y = altitude, col = factor(sta_id))) +
  # geom_vline(data = twl, aes(xintercept = twilight, linetype = ifelse(rise, "dashed", "solid"), color = "grey"), lwd = 0.1) +
  theme_bw() +
  scale_colour_manual(values = col) +
  scale_y_continuous(name = "Altitude (m)")

ggplotly(p, dynamicTicks = T) %>% layout(showlegend = F)

Wind assistance

Show code
  fun_marker_color <- function(norm){
    if (norm < 20){
      "darkpurple"
    } else if (norm < 35){
      "darkblue"
    } else if (norm < 50){
      "lightblue"
    } else if (norm < 60){
      "lightgreen"
    } else if (norm < 80){
      "yellow"
    } else if (norm < 100){
      "lightred"
    } else {
      "darkred"
    }
  }
  fun_NSEW <- function(angle){
    angle <- angle  %% (pi* 2)
    angle <- angle*180/pi
    if (angle < 45/2){
      "E"
    } else if (angle < 45*3/2){
      "NE"
    } else if (angle < 45*5/2){
      "N"
    } else if (angle < 45*7/2){
      "NW"
    } else if (angle < 45*9/2){
      "W"
    } else if (angle < 45*11/2){
      "SW"
    } else if (angle < 45*13/2){
      "S"
    }else if (angle < 45*15/2){
      "SE"
    } else {
      "E"
    }
  }

  sta_duration <- unlist(lapply(static_prob_marginal,function(x){as.numeric(difftime(metadata(x)$temporal_extent[2],metadata(x)$temporal_extent[1],units="days"))}))

  m <-leaflet(width = "100%") %>%
    addProviderTiles(providers$Stamen.TerrainBackground) %>%  addFullscreenControl() %>%
    addPolylines(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#808080", weight = 3) %>%
    addCircles(lng = shortest_path$lon, lat = shortest_path$lat, opacity = 1, color = "#000", weight = sta_duration^(0.3)*10)

  for (i_s in seq_len(grl$sz[3]-1)){
    if (grl$flight_duration[i_s]>5){
      edge <- which(grl$s == shortest_path$id[i_s] & grl$t == shortest_path$id[i_s+1])

      label = paste0( i_s,': ', grl$flight[[i_s]]$start, " - ", grl$flight[[i_s]]$end, "<br>",
                      "F. dur.: ", round(grl$flight_duration[i_s]), ' h <br>',
                      "GS: ", round(abs(grl$gs[edge])), ' km/h, ',fun_NSEW(Arg(grl$gs[edge])),'<br>',
                      "WS: ", round(abs(grl$ws[edge])), ' km/h, ',fun_NSEW(Arg(grl$ws[edge])),'<br>',
                      "AS: ", round(abs(grl$as[edge])), ' km/h, ',fun_NSEW(Arg(grl$as[edge])),'<br>')

      iconArrow <- makeAwesomeIcon(icon = "arrow-up",
                                   library = "fa",
                                   iconColor = "#FFF",
                                   iconRotate = (90 - Arg(grl$ws[edge])/pi*180) %% 360,
                                   squareMarker = TRUE,
                                   markerColor = fun_marker_color(abs(grl$ws[edge])))

      m <- m %>% addAwesomeMarkers(lng = (shortest_path$lon[i_s] + shortest_path$lon[i_s+1])/2,
                                   lat = (shortest_path$lat[i_s] + shortest_path$lat[i_s+1])/2,
                                   icon = iconArrow, popup = label)
    }
  }
  m

Histogram of Speed

Show code
edge <- t(graph_path2edge(path_sim$id, grl))
nj <- ncol(edge)
nsta <- ncol(path_sim$lon)

speed_df <- data.frame(
  as = abs(grl$as[edge]),
  gs = abs(grl$gs[edge]),
  ws = abs(grl$ws[edge]),
  sta_id_s = rep(head(grl$sta_id,-1), nj),
  sta_id_t = rep(tail(grl$sta_id,-1), nj),
  flight_duration = rep(head(grl$flight_duration,-1), nj),
  dist = geosphere::distGeo(
    cbind(as.vector(t(path_sim$lon[,1:nsta-1])), as.vector(t(path_sim$lat[,1:nsta-1]))),
    cbind(as.vector(t(path_sim$lon[,2:nsta])),   as.vector(t(path_sim$lat[,2:nsta])))
  ) / 1000
) %>% mutate(
  name = paste(sta_id_s,sta_id_t, sep="-")
)

plot1 <- ggplot(speed_df, aes(reorder(name, sta_id_s), gs)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot2 <- ggplot(speed_df, aes(reorder(name, sta_id_s), ws)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot3 <- ggplot(speed_df, aes(reorder(name, sta_id_s), as)) + geom_boxplot() + theme_bw() +scale_x_discrete(name = "")
plot4 <- ggplot(speed_df, aes(reorder(name, sta_id_s), flight_duration)) + geom_point() + theme_bw() +scale_x_discrete(name = "")

subplot(ggplotly(plot1), ggplotly(plot2), ggplotly(plot3), ggplotly(plot4), nrows=4, titleY=T)